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ABSTRACT 


Manned  aircraft  that  are  intended  for  surveillance  or 
to  complete  a  bombing  mission  will  very  likely  be  engaged 
by  surface  to-air-missiles  having  guidance  systems  based  on 
infrared  (IR)  technology.  The  objective  of  this  study  was  to 
characterize  via  simulation  the  amount  of  "cover"  that  can 
be  obtained  by  dropping  from  a  pre- launched,  unmanned 
tactical  air  launched  decoy  (TALD)  a  sequence  of  pyrophoric 
materials  to  create  an  IR  cloud,  analogous  to  the 
interference  created  by  microwave  chaff,  that  would  protect 
the  manned  aircraft  from  the  missile.  The  performance 
analysis  is  based  on  a  simple  reticle  based  model  in  which 
the  two-dimensional  (2D)  image  is  reduced  to  either  a 
composite  signal,  created  by  the  aircraft,  or  a  composite 
noise,  created  by  the  pyrophoric  expandable.  The  analysis 
leads  to  a  computer  simulation  model  producing  time  and 
space  dependent  signal-to-noise  ratios.  It  is  demonstrated 
that  the  simulation  model  can  answer  questions  such  as  how 
long  the  materials  need  to  burn,  how  much  intensity  is 
needed,  what  wavelength  range  is  most  effective,  which 
pyrophoric  packets  should  be  dropped,  and  how  many.  A 
visual  model  of  the  time  dependent  IR  pyrophoric  cloud  has 
also  been  created. 
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I 


INTRODUCTION 


A.  BACKGROUND 

The  most  effective  anti-aircraft  weapon  ever  created 
to-date  is  the  infrared  (IR)  guided  missile.  Since  their 
introduction  into  operational  use  in  the  early  1950s,  IR 
missiles  have  far  exceeded  the  expectations  their  designers. 
The  enabling  technology  for  practical  IR  missiles  came  from 
WWII  research  on  detector  materials  that  could  be  used  to 
fashion  detectors  sensitive  in  the  IR  bands.  These  initial 
materials  were  sensitive  to  radiation  in  the  near-IR  (1-2 
micron  wavelength)  band.  Over  the  last  20  years,  advanced 
designs  have  taken  advantage  of  better  detector  materials, 
moving  the  bandpass  from  the  near-IR  into  the  mid-IR  (3-5 
micron  band)  ,  where  engine  plxime  provides  a  superior  signal 
for  all-aspect  engagement.  The  three  to  five  micron  band 
also  has  less  atmospheric  attenuation  and  clutter  factors. 

The  need  to  protect  aircraft  from  attack  by  effective, 
easily  launched  IR  homing  missiles  led  to  the  development  of 
pyrotechnic  flares.  The  flare  needs  to  have  a  jamming-to- 
signal  (J/S)  ratio  of  greater  than  one  to  one  versus  the 
protected  aircraft's  signature  to  lure  a  missile  away  from 
the  aircraft .  Since  the  observed  surface  area  of  the  IR 
signature  of  the  flare  is  quite  small,  a  pyrotechnic  flare 
has  to  operate  at  very  high  temperatures  to  match  or  exceed 
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the  in-band  aircraft  signature.  But  high  temperatures  shift 
the  wavelength,  and  we  are  forced  by  the  laws  of  physics  to 
increase  the  size  of  the  flare  to  provide  enough  energy  at 
the  correct  wavelength.  The  flare's  intensity  is  presumed  to 
rapidly  rise  to  several  orders  of  magnitude  higher  than  that 
of  the  target  signature  and  ignite  near  the  target  aircraft, 
inside  the  field  of  view  (FOV)  of  missile's  seeker. 

B.  APPROACH 

In  this  thesis  the  amount  of  protection  that  can  be 
obtained  by  dropping  a  sequence  of  pyrophoric  materials  to 
create  an  IR  chaff  cloud  will  be  simulated.  A  simple 
scenario  will  be  defined  in  which  the  packets  are  dropped 
from  a  pre-launched  unmanned  tactical  air  launched  decoy 
(TALD)  ,  and  the  manned  aircraft  follows,  at  a  later  time, 
close  to  the  same  path  as  the  decoy,  taking  advantage  of  the 
interfering  effect  created  by  the  pyrophoric  clouds. 

The  heat-seeking  missile  is  .  assumed  to  be  in  a 
prelaunch  status  so  that  missile  aero-dynamics  are  not  part 
of  this  initial  study.  The  imaging  capability  of  the  missile 
is  assumed  to  be  defined  in  teirms  of  a  definable  number  of 
pixels  within  the  known  field-of-view.  Also  assumed  known  is 
the  detector  sensitivity  and  detectivity  over  one  or  more 
bands  of  wavelengths.  The  image  of  the  aircraft  is  defined 
in  terms  of  a  distribution  of  hot  spots.  Any  interference 
provided  by  the  IR  signature  of  the  TALD  aircraft  are 
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ignored.  A  visual  model  for  expenditure  of  the  IR  chaff 
cloud  is  created.  The  sequence  and  timing  for  the  dropping 
of  the  packets  are  assumed  to  be  known.  Also,  both  the  speed 
of  the  manned  aircraft  and  the  decoy  are  assumed  known. 

The  analysis  of  the  problem  can  answer  questions  such 
as  how  long  the  materials  need  to  burn,  how  much  intensity 
is  needed,  what,  wavelength  range  and  distribution  are  most 
effective,  the  rate  at  which  pyrophoric  packets  should  be 
dropped  and  how  many,  and  lastly,  quantitatively  what  kind 
of  signal-to-noise  ratio  (S/N)  deterioration  can  be 
extracted.  The  mathematical  software  MATLAB  will  be  used  to 
simulate  the  pyrophoric  material,  provide  image  snapshots  in 
time  of  the  pyrophoric  material,  and  calculate  the  resulting 
S/N. 

Some  basic  IR  theoiry  such  as  Plank's  Law,  Wien's  Law, 
and  the  Stefan-Boltzmann  Law,  are  reviewed  in  Chapter  II. 
Radiometric  definitions,  atmospheric  absorption,  and  general 
characteristics  of  the  different  types  of  IR  detectors  are 
also  discussed  in  Chapter  II.  The  design  of  the  pyrophoric 
model  is  discussed  in  Chapter  III.  The  design  of  the 
aircraft  plume  model  is  addressed  in  Chapter  IV.  The 
integration  of  the  pyrophoric  and  pl\ime  model  is  discussed 
in  Chapter  V.  The  first  part  of  the  simulation  code  and  some 
S/N  curves  are  presented  in  Chapter  VI.  The  second  part  of 
the  simulation  code  and  the  pyrophoric  and  pl\ame  images  are 
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presented  in  Chapter  VII .  The  use  of  multiple  pyrophoric 
expendables  is  discussed  in  Chapter  VIII.  Finally, 
conclusions  from  the  study  of  the  pyrophoric  model  and 
future  enhancements  in  modeling  are  presented  in  Chapter  IX. 
Also,  in  Appendix  A,  B  and  C  is  a  brief  description  of  the 
symbols  and  the  quantities  that  are  used  during  the  analysis 
or  in  the  simulation  code,  the  MATLAB  codes,  and  the 
mathematical  derivations,  respectively. 
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II.  BACKGROUND 


A.  INTRODUCTION 

In  the  year  1800,  Sir  William  Herschel,  the  royal 
astronomer  to  the  King  of  England,  was  conducting  an 
experiment  with  a  prism  in  sunlight.  The  prism  spread  the 
sun's  rays  into  a  spectrum  from  violet  to  red.  Herschel 
placed  a  thermometer  in  the  violet  color  and  recorded  the 
temperature.  He  moved  the  thermometer  through  the  colors 
from  blue  to  red  and  noticed  that  the  temperature  increased 
progressively.  He  then  moved  the  thermometer  beyond  the  red 
end  of  the  visible  region  and  the  temperature  continued  to 
increase.  Thus,  he  found  energy  beyond  the  red;  this  energy 
has  come  to  be  known  as  infrared  (IR) . 

A  portion  of  the  electromagnetic  spectrum  is  indicated 
in  Figure  2.1,  All  electromagnetic  radiation  obeys  similar 
laws  of  reflection,  refraction,  diffraction,  and 
polarization.  The  velocity  of  propagation  is  the  same  for 
all.  They  differ  from  one  another  only  in  wavelength  and 
frequency.  The  portion  of  the  spectrxim  that  includes  the 
infrared  is  depicted  in  greater  detail  in  the  lower  part  of 
Figure  2.1,  where  the  infrared  region  is  bounded  on  the 
short-wavelength  side  by  visible  light  and  on  the  long- 
wavelength  side  by  microwaves.  It  is  convenient  to  subdivide 
the  infrared  region  into  several  parts .  There  are  no  exact 
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designations  for  the  separation  of  infrared  bands.  The  most 
universally  accepted  designations  today  are  as  follows: 
visible  light  up  to  3  urn  is  designated  near  infrared  (NIR) , 
inf^si^sd  (MIR)  is  from  3  to  6  um,  and  far  infrared 
is  from  6  to  15  um.  Beyond  15  um  is  the  extreme 
infrared  (XIR) . 


j  visaiit 
!  V  t  CYOU 


Near  infrared 


xHjz 


0.6  QJ6  1 
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Middle  infrared  [  Far  infrared  [  Extreme  infrared  ^ 

- "i  I  j  I  I  1 1  I  < 

^  *  6  8  10  B  a)  30 


Figure  2.1  The  Electromagnetic  Spectrum.  (From  Ref  [1] ) 


B .  RADIOMETRIC  SYMBOLS 

All  objects  with  nonzero  temperature  will  emit 
electromagnetic  radiation.  Some  objects  prove  to  be  better 
radiators  than  others.  The  (hypothetical)  best  radiator  is 
one  that  obeys  a  set  of  classical  equations  called  the 
blackbody  radiation  equations .  Some  of  te2rms  used  to 
describe  IR  sources  and  radiation  are  defined  in  Table  2.1. 
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Table  2.1  Symbols,  Descriptions,  and  Units 


Syntool 

'  Term  .  ■  /  ' 

Description 

Unit 

U 

Radiant 

energy 

Energy  transferred  by 

electromagnetic  waves 

Joule 

P 

Radiant  flux 

Rate  of  transfer  of 

radiant  energy 

Watt  (W) 

W 

Radiant 

emittance 

Radiant  fliix  emitted  per 

unit  area  of  a  source 

W  •  cnf^ 

J 

Radiant 

intensity 

Radiant  flux  per  unit 

solid  angle 

W-cm^ 

N 

Radiance 

Radiant  flux  per  unit 

solid  angle  per  unit  area 

W  •  cm 

H 

Irradiance 

Radiant  flux  incident  per 

unit  area 

Wx 

Spectral 

radiant 

emittance 

Radiant  emittance  per  unit 

wavelength  interval  at  a 

particular  wavelength 

W  -cm 

E 

Emissivity 

Ratio  of  radiant  emittance 

of  a  source  to  that  of  a 

blackbody  at  the  same 

temperature 

Dimensionless 
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C.  BLACKBODY  RADIATORS 

The  spectral  density  distribution  of  the  power  radiated 
into  a  hemisphere  from  a  blackbody  radiator  per  unit  area  of 
the  source  is  given  by  Planck's  Law  [2]: 


iTlhc^ 


1 


where 

-  Wx  is  the  spectral  radiant  emittance  of  the  source  with 

units  of  W-cm'^-(im‘^ 

-  h  is  Planck's  constant  (h=6 . 625xl0'^‘  joules-s) 

-  c  is  the  vacuum  velocity  of  light  {c=3xl0*  m-s'^) 

-A,  is  the  wavelength  of  the  radiation  {|am  or  nm) 

-  k  is  Boltzman's  constant  (k=1.38xl0"”  joules-K'^) 

-  T  is  the  temperature  of  the  source  in  K 

-  c,=  27rhc'=3.74xl0'  W-cm'*-Mm' 

-  Cj=  hc/k=l.  439x10"  Jtm-K 


Plots  of  the  spectral  radiant  emittance  for  a  typical  range 
of  temperatures  encountered  in  our  environment  are  shown  in 
Figure  2.2.  The  spectral  radiant  emittance  is  the  spectral 
power  distribution  normalized  by  the  area  of  the  source. 
Note  that  most  of  the  energy  is  in  the  infrared  portion  of 
the  spectrum.  As  the  temperature  of  the  source  increases,  we 
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/ 

note  that  the  area  under  the  curve  increases,  and  the 
location  of  the  peak  shifts  toward  the  shorter  wavelengths. 


Figure  2.2  Spectral  Radiant  Emittance  vs.  Wavelength  for 
Blackbodies  of  Various  Tenperatures .  (From  Ref  [2]) 


The  location  of  the  peak  of  the  spectral  radiant  emittance 
is  foiind  by  taking  the  derivative  of  Planck's  Law  with 
respect  to  wavelength.  It  can  be  shown  that  the  result  gives 
Wien's  Law,  the  wavelength  of  the  peak  spectral  radiant 
emittance  : 


2.896x10^ 

T[K] 


[jum] 


{II -2) 
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Knowing  the  temperature  of  the  source,  we  can  calculate  the 
location  of  the  peak  emittance  and,  as  shown  in  Figure  2.2, 
higher  temperatures  imply  lower 


The  value  of  the  peak  spectral  radiant  emittance  is  found  by 
sxibstituting  (II-2)  into  (II-l)  to  obtain 


=  ^bT= 


where 


b  =  1.286 X 10"*^ W  •  cm-^  • 


The  total  emittance  W  is  the  total  power  emitted  into  a 
hemisphere  at  all  wavelengths  per  unit  area  of  the  source. 
It  is  found  by  integrating  Planck's  Law  over  all 
wavelengths , 

OO 

(77-4) 

0 

The  result  of  this  integration  is 


W  = 


ISc^h^ 


=  cfr* 


where 


a  =  5.67xlO"‘^W 


(fl-5) 


which  represents  the  Stefan-Boltzmann  Law  and  shows  that  the 
total  power  per  unit  area  of  the  source  increases  as  the 
fourth  power  of  the  temperature. 
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We  frequently  pass  electromagnetic  radiation  through  a 
spectral  filter  (like  the  atmosphere)  with  a  spectral 
response  of  the  form  S{X)  .  The  emittance  that  is  passed  by 


such  a  filter  is 

w(;ii  ,^2 ) = J  s{Xj^;^dX  {n  -  6) 

Figure  2.3  is  an  illustration  of  one  such  graph  where  the 
vertical  axis  is  used  to  compute  the  fraction  of  the  total 
emittance  that  lies  between  a  wavelength  of  0  |tm  and  some 
upper  wavelength  divided  by  the  total  emittance  over  all 
wavelengths . 


Figure  2.3  W(0,X,)  Relative  to  Total  Emittance.  (Ref  [2]) 


D.  RADIATION  FROM  OBJECTS 

In  order  to  calculate  the  radiation  emitted  from  a 
source,  we  model  the  source  as  either  a  uniformly  radiating 
spherical  source  (e.g,  the  sun  or  a  light  bulb)  or  as  a 
planar  source (e.g,  the  exhaust  port  of  a  jet  engine) . 

1.  Uniform  Spherical  Source 

For  a  spherical  source  with  a  known  temperature ( in  K) , 
an  emitting  area  A,  and  an  emissivity  (usually  assximed  to  be 
less  than  one;i.e.,  a  graybody  object  [2]),  we  can 
calculate  the  radiant  emittance  of  the  source  as 

w  =  s&r''  (//-7) 


Then  we  can  find  the  total  power  from  the  source  as 


P  =  W-A  (/7_8) 

and  the  radiant  intensity  by  realizing  that  it  is  emitting 

uniformly  into  a  sphere  of  An  steradians: 


(II -9) 


2 .  Planar  Source 

For  a  planar  source,  we  cannot  assume  that  the 
radiation  is  symmetrical  so  the  power  is  radiated 
nonuni formly  from  such  a  source.  A  standard  procedure  is  to 
model  the  source  as  a  Lambertian  emitter,  where  both  power 
and  radiant  intensity  will  vary  with  the  observation  angle. 
For  such  a  source  it  can  be  shown  that  the  radiance  is  [2] 

N  =  =  (//-lO) 

m  71 
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E.  TRANSMISSION  OF  IR  THROUGH  THE  EARTH'S  ATMOSPHERE 

Most  infrared  systems  must  view  their  targets  through 
the  earth's  atmosphere.  Before  it  reaches  the  infrared 
sensor,  the  radiant  flux  from  the  target  is  selectively 
absorbed  by  several  of  the  atmospheric  gases,  and  scattered 
away  from  the  line  of  sight  by  small  particles  suspended  in 
the  atmosphere. 

When  the  particles  are  small  compared  with  the 
wavelength  of  the  radiation,  the  process  is  known  as 
Rayleigh  scattering  and  exhibits  a  dependence.  For  larger 

particles,  the  scattering  is  independent  of  wavelength. 
Scattering  by  gas  molecules  in  the  atmosphere  is,  therefore, 
negligibly  small  for  wavelengths  longer  than  2  -urn.  Smoke  and 
light  mist  particles  are  also  usually  small  with  respect  to 
infrared  wavelengths,  and  infrared  radiation  can,  therefore, 
penetrate  further  through  smoke  and  mists  than  visible 
radiation.  However,  rain,  fog  particles,  and  aerosols  are 
larger  and,  consequently,  scatter  infrared  and  visible 
radiation  to  a  similar  degree. 

In  the  infrared  portion  of  the  spectriim,  the  absorption 
process  poses  a  far  more  serious  problem  than  does  the 
scattering  process.  The  spectral  transmittance  measured  over 
a  6000  ft  horizontal  path  at  sea  level  is  shown  in  Figure 
2.4.  The  molecule  responsible  for  each  absorption  band, 
either  water  vapor,  carbon  dioxide,  or  ozone,  is  shown  in 
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the  lower  part  of  the  figure.  Inspection  of  Figure  2,4 
reveals  the  presence  of  atmospheric  windows,  i.e,  regions  of 
reduced  atmospheric  attenuation. 

IR  detection  systems  are  designed  to  operate  in  these 
windows.  Combinations  of  detectors  and  spectral  bandpass 
filters  are  selected  to  define  the  operating  region  to 
conform  a  window  to  maximize  performance  and  minimize 
background  contributions. 


twit  1 1  It'  t  t  It  t 

H,0  cc4^  “b,  “2°  Hrf)  CO,  CO2 


ABSORBING  MOLECULE 


Figure  2.4  Transmittance  of  Atmosphere  over  one  Nautical 
Mile  (NM)  Sea  Level  Path  (Infrared  Region).  (From  Ref  [3]) 
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F.  IR  DETECTORS 

The  detector  is  the  heart  of  every  IR  system  because  it 
converts  scene  radiation  into  some  other  measurable  form; 
this  can  be  an  electrical  current,  or  a  change  in  some 
physical  property  of  the  detector. 

In  general,  an  infrared  detector  in  a  particular  set  of 
operating  conditions  is  characterized  by  two  performance 
measures:  the  responsivity  R  and  the  specific  detectivity 
D* .  The  responsivity  is  the  gain  of  the  detector  expressed 
in  volts  of  output  signal  per  watt  of  input  signal.  The 
specific  detectivity  is  the  detector  output  signal-to-noise 
ratio  for  one  watt  of  input  signal,  normalized  to  a  unit 
sensitive  detector  area  and  a  vinit  electrical  bandwidth.  The 
different  types  of  IR  detectors  may  be  divided  into  two 
broad  classes,  namely  thermal  detectors  and  photon,  or 
quant\im,  detectors. 

1 .  Thermal  Detectors 

Because  thermal  detectors  have  been  used  since 
Herschel ' s  discovery  of  the  infrared  portion  of  the 
spectrum,  it  is  appropriate  to  consider  them  first.  They  are 
distinguished  as  a  class  by  the  observation  that  the  heating 
effect  of  the  incident  radiation  causes  a  change  in  some 
physical  property  of  the  detector. 

Since  most  thermal  detectors  do  not  require  cooling, 
they  have  found  almost  universal  acceptance  in  certain  field 
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applications  in  which  it  is  impractical  to  provide  such 
cooling.  Because  (theoretically)  they  respond  equally  to  all 
wavelengths,  thermal  detectors  are  often  used  in 
radiometers .  The  time  constant  of  a  thermal  detector  is 
usually  a  few  milliseconds  or  longer,  so  that  they  are 
rarely  used  in  search  systems  or  in  any  other  application  in 
which  high  data  rates  are  required. 

2 .  Photon  Detectors 

Most  photon  detectors  have  a  detectivity  that  is  one  or 
two  orders  of  magnitude  greater  than  that  of  thermal 
detectors.  This  higher  detectivity  does  not  come  for  free, 
however,  since  many  photon  detectors  will  not  function 
unless  they  are  cooled  to  cryogenic  temperatures .  Because  of 
the  direct  interaction  between  the  incident  photons  and  the 
electrons  of  the  detector  material,  the  response  time  of 
photon  detectors  is  very  short;  most  have  time  constants  of 
a  few  microseconds  rather  than  the  few  milliseconds  typical 
of  theional  detectors. 

Finally,  the  spectral  response  of  photon  detectors, 
unlike  that  of  thermal  detectors,  varies  with  wavelength. 
Curves  of  D*  versus  wavelength  are  shown  in  Figure  2 . 5  for 
different  types  of  photon  detectors. 
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Figure  2.5  D*  vs.  X  for  Representative  Detectors  (From  Ref 
[2]) 
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III.  PYROPHORIC  MODEL 

A.  GENERAL  DESIGN  REQUIREMENTS 

Infrared  decoys  are  used  in  modern  times  to  protect 
combat  systems  from  IR  tracking  threats .  The  most  common 
examples  of  IR  decoys  are  IR  flares .  IR  decoys  are  used  to 
protect  aircraft  from  heat-seeking  missiles.  Some  key 
requirements  for  the  design  specifications  of  an  expendable 
decoy  are  discussed  in  the  following. 

1.  Peak  Intensity 

Peak  intensity  is  normally  the  most  important 
requirement.  IR  decoys  must  radiate  with  sufficient 
intensity  and  at  least  exceed  the  intended  target's  radiant 
intensity  in  the  band  of  interest.  This  can  be  accomplished 
by  controlling  the  decoy  temperature  or  by  the  use  of 
selective  emitters  that  have  a  higher  emissivity  in  the  band 
of  interest.  A  common  design  problem  is  illustrated  in  the 
spectral  intensity  versus  wavelength  plot  of  Figure  3.1.  The 
relative  spectral  distribution  between  a  decoy  (a  small  hot 
source)  and  an  aircraft  target  (a  large,  relatively  cool 
source)  is  shown.  From  the  viewpoint  of  the  decoy  designer, 
it  is  advantageous  if  a  significantly  higher  percentage  of 
the  optical  signal  at  the  detector,  within  the  '  decoy  band' 
is  due  to  the  decoy  rather  than  the  target.  Finally,  the 
peak  intensity  is  the  primary  driver  for  the  decoy  weight. 
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volume,  and  cost  that  must  be  smaller,  lighter,  and 
cheaper, respectively,  than  the  protected  target. 


2.  Flare  Intensity  Rise  Time 

A  decoy  must  persist  long  enough  to  ensure  that  no 
possibility  of  target  reacquisition  remains.  Consequently, 
it  must  maintain  a  credible  signature  until  the  original 
target  is  no  longer  in  the  threat  f ield-of-view.  If  this  is 
not  the  case,  it  is  necessary  to  deploy  a  second  flare. 
Additionally,  the  decoy  must  achieve  an  effective  intensity 
quickly  enough  to  capture  the  seeker  before  leaving  the 
threat  field-of-view.  The  diameter  of  the  threat  field-of- 
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view  at  the  time  of  decoy  deployment  is  usually  less  than 
200  m.  This  means  that  effective  operational  intensity 
levels  must  be  achieved  in  a  fraction  of  a  second  [Ref  4]  . 
The  exact  value  of  the  time  to  rise  to  peak  radiant 
intensity  depends  of  the  chemical  composition  and  packing. 
Exact  values  are  difficult  to  find  in  the  open  literature. 

B.  DEVELOPMENT  OF  THE  PYROPHORIC  MODEL 

The  pyrophoric  model  must  simulate  changes  in  both  the 
time  and  space  domains .  The  modeling  procedure  for 
calculating  the  in-band  radiant  intensity  is  summarized  in 
Figure  3.2. 


Figure  3.2  Calculation  of  In-Band  Radiant  Intensity. 
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1.  Time  Domain  Variations 

The  emittance  over  all  wavelengths  for  any  hot  object 
follows  Planks  radiation  equation  is  given  from  the  Stefan- 
Boltzmann  Law 


W{A.T)„=scff,/ 

where  the  Stefan  Boltzmann  radiation  constant  is  [Ref  2] 

(j  =  5.67  X 10"^^ W  •  cm-^  °K-^ 

e  is  the  emissivity  of  the  pyrophoric  and  Tp^  a  specified 
temperature.  The  approximation  of  a  graybody  model  is  used 
to  approximate  the  in-band  radiant  emittance  for  these 
expendables  [Ref  4] . 

The  M-file  plankpf.m  recorded  for  reference  in  Appendix 
B  is  used  to  calculate  for  a  particular  temperature  the 
integrated  emittance  over  a  band  of  wavelengths.  The  3-to  5- 
micron  band  is  very  common  for  a  detector  because  this  band 
has  less  atmospheric  attenuation.  The  emittance  W(X^,Xj)pj 
that  lies  inside  the  detector's  band  is  obtained  using  the 
M-file  program  cited  above.  The  fraction  of  the  emitted 
power  that  lies  between  two  wavelengths  and  X^  is 


The  pyrophoric  source  is  modeled  as  a  uniformly 
radiating  spherical  source  with  radius  r^  at  the  instant 
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for  which  a  maximum  occurs  in  the  radiant  intensity  t^.  The 
associated  area  is 


Apj  =47rrt 


2 

MAX 


(III -3) 


The  peak  radiant  intensity  of  an  isotropic  radiator  is 
[Ref  2] 


(m-A) 


The  total  radiant  intensity  that  lies  inside  the  detector's 
band  is 

~  J MAX  ■  Vpf  ~ 


From  the  input  t^  for  the  pyrophoric  burn,  the  time 
radiant  intensity  distribution  of  the  pyrophoric  burn  can  be 
predicted.  Two  operational  properties  are  desirable. 
Ideally,  the  pyrophoric  time  function  will  achieve  peak 
intensity  quick  enough  to  capture  missile  seeker  attention 
before  leaving  the  f ield-of-view.  Furthermore,  it  should 
persist  long  enough  to  ensure  that  no  possibility  of  target 
reacquisition  remains.  The  probability  density  function 
(pdf)  of  the  standard  Gamma  distribution  [Ref  5]  has 
sufficient  flexibility  to  model  these  two  basic  properties. 
The  flow  chart  that  s\ammarizes  the  development  of  the 
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temporal  variation  in  radiant  intensity  distribution  of  the 
pyrophoric  burn  is  shown  in  Figure  3.3. 


/Set 

Pyrophoric 

Inputs 
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C  Create  the  time  radiant  intensity 
distribution  of  pyrophoric  bum 


3.3  The  Time  Radiant  Intensity 
Structure  of  Pyrophoric  Burn. 


Distribution 
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other  inputs  such  as  that  corresponding  to  initial  time  t^ 
and  final  time  t^  for  the  simulation  to  be  performed  will 
lead  to  the  development  of  a  time-function-array  t. 

The  Gamma  distribution,  which  is  a  continuous  random 
variable  (RV) ,  is  defined  as 


t>0 

a>0  ,  /3>0 


(III -6) 


Because  of  interest  in  modeling  the  teir^oral  variation  in 
radiant  intensity,  the  t  parameter  is  used  to  define  time 
dependence.  The  standard  Gamma  distribution  [Ref  5]  has  P=1 
so  the  pdf  of  a  RV  modeled  as  a  standard  Gamma  RV  is 


r>0 

a>0 

In  Appendix  C,  a  short  proof  is  presented  that  demonstrates 
that  the  peak  in  the  distribution  satisfies  the  condition 

(ffi-8) 

From  this  relation  the  values  for  the  parameter  a  for 
different  values  of  the  pyrophoric  peak  times  are 
predicted.  Using  the  M-file  gaunmaz.m,  recorded  for  reference 
in  Appendix  B,  we  calculate  the  Gamma  function  r(a)  and 
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Figure  3.4  The  Time  Radiant  Intensity  Distribution  of 
Pyrophoric  Burn 
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2 .  Space  Domain 

The  flow  chart  that  is  followed  to  develop  the 
radiant  intensity  distribution  of  pyrophoric  burn  in  the 
space  domain  is  shown  in  Figure  3.5. 


Figure  3.5  Normalized  Space  Distribution  Structure 
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It  is  assTjmed  that  the  spatial  distribution  for  each 
instant  of  the  pyrophoric  burn  is  a  Gaussian  distribution. 
Generally,  a  continuous  random  variable  X  is  said  to  have  a 


Gaussian  distribution  if  the  pdf  is 


f{x\n,a)  = 


‘^2-71  • 


2-cr^ 


-oo<r<oo 


{iii-m) 


where  |X  is  the  mean  and  a  is  the  standard  deviation.  One 
Gaussian  distribution  is  shown  in  Figure  3.6.  The  Gaussian 
pdf  is  symmetric  about  the  mean  and  the  standard  deviation 
is  the  distance  from  the  mean  to  the  point  at  which  the 
slope  of  the  curve  changes. 


Figure  3.6  Graph  of  a  Gaussian  Distribution 


In  Appendix  C  a  proof  is  presented  that  demonstrates  that 
for  a  Gaussian  random  variable  Z,  the  area  under  the  curve 
from  zero  to  infinity  is  given  by  equation  (C-12),  which  is 
reproduced  here  for  convenience: 


(ffl-ii) 


From  the  inputs  t^,  r^,  and  the  pyrophoric  time  function 

array  t,  the  following  parameters  can  be  obtained: 


(cons tmt  pyrophoric  velocity) 


{temporal  dependent  pyrophoric  radius) 


^NORM  =  {normalized  pyrophoric  radius) 

^MAX 


(///-12) 

(///-13) 


(/7/-14) 


At  each  instant  of  the  pyrophoric  burn,  it  is  ass-umed  that 
the  spatial  radiant  intensity  function  is  a  Gaussian 
distribution  with  a  mean  value 


(///-15) 

From  (III-ll)  ,  the  area  under  each  curve  from  zero  to 
infinity  is 


, ,  _  ^MAXt  _  ^  _  tj 


•MAX 


U  't 

P 
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(///-16) 


which  is  positive  and  less  than  one. 

For  each  element  of  the  array  t,  the  value  of  the 
temporal  distribution  is  given  by 


(lll-ll) 


Now  in  equation  (III-IO)  we  use  x=  r^,^  and  (A=  t/tp  to  obtain 

the  normalized  space  distribution  in  the  pyrophoric  burn 
as 


fin 


NORM’ 


{ill -is) 


The  next  step  is  to  use  1/Q(t)  to  normalize  the  corresponding 
equation  (III-18) .  The  result,  which  is  a  pdf  defined  over 
positive  r^oRM'  an  area  under  the  curve  equal  to  one. 
This  result  is  multiplied  by  equation  (III-9)  to  obtain 


'NORM’  /f 


-fin 


NORM’* 


{in -19) 


where  the  space  integration  of  the  space-time  distribution 
satisfies  the  relation 
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*NORM^  /f 


\fr  =  7(r) 


V 


This  forces  the  resulting  space-time  distribution  j  to  have 
a  space  integrated  area  equal  to  J{t) .  One  example  is  shown 
in  Figure  3.7  where  successive  curves  represent  time  steps 
of  0.2s. 

s 

x10 


Figure  3.7  Pyrophoric  Space  Distribution 
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The  last  major  step,  see  the  flowchart  in  Figure  3.5,  is  to 
create  the  space  domain  pyrophoric  integrated  view-model. 
The  general  idea  is  to  pick  a  curve  from  Figure  3.7  which 
corresponds  to  one  instant  of  time  t,  as  illustrated  in 
Figure  3.8. 


^Nom 

3.8  Pyrophoric  Space  distribution  for  one  instant  of 

time  t . 

As  represented  on  Figure  3.9,  the  pyrophoric  source  is 
modeled  as  spherical.  It  is  assumed  that  an  external  remote 
observer  will  'view'  the  source  and  perceive  radiation 
integrated  along  a  thin  line  x.  The  total  radiant  intensity 
distribution  at  a  specific  instant  of  time  will  consist  of 
the  integration  of  radiant  intensity  over  each  distance  |  rj, 
where  |  r.|  is  the  distance  between  the  center  of  the 
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distribution  and  some  point  along  the  line  x.  The  line  x  is 
shifted  by  |  rj  from  the  distribution  center. 


Figure  3.9  The  integrated  pyrophoric  model  projection 


Analytically,  for  each  distance  from  the  center  of  the 
distribution  |  rJ,  a  fixed  number  of  uniformly  spaced  points 
are  selected  along  the  line  x  from  x=0  to  x=x^.  The 
distance  |  rJ  that  connects  each  of  these  points  with  the 
center  of  the  distribution  is  given  by 


where  |  x^|  the  distance  along  the  line  x.  In  the  computer 
code,  an  index  k  is  generated  as 


(///-2l) 


so  the  integrated-view  radiant  intensity  for  the  specific 
I  rj,  i.e,  one  point  is  given  by  jt(k).  The  siainmation  of  all 
the  radiant  intensity  values  for  each  |  r^|  give  the 
integrated-view  radiant  intensity 


^  [in -22) 

which  is  easily  computed.  Alternatively,  an  equivalent  but 
more  analytical  representation  is 


{in -23) 


Typical  distributions  of  integrated-view  radiant  intensities 
for  pyrophoric  burns  lasting  eight  seconds  are  shown  in 
Figure  3.10.  The  horizontal  axis  of  Figure  3.10  is  specified 
in  terms  of  the  index  variable  of  the  array  r„.  It  is 

observed  that  the  curve  with  the  maximum  area  corresponds  to 
the  instant  t=t  . 

P 
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IV.  PLtmE  MODEL 


A.  GENERAL  DESIGN  REQUIREMENTS 

The  pliime  of  the  aircraft  can  dominate  the  IR 
signature.  This  is  no  doubt  due  to  the  considerable  radiant 
energy  created  from  the  combustion  process .  Measurements  of 
the  radiation  from  military  a  aircraft's  pliime  are  highly 
classified.  Fortunately,  it  is  not  too  difficult  to  make 
order-of-magnitude  calculations  of  the  radiation  by  using 
only  temperature  and  dimensional  information  that  can  be 
found  in  the  open  literature  [Ref  4] . 

During  afterburning,  which  results  in  an  increase  in 
thrust,  the  rate  at  which  fuel  is  consumed  increases 
drastically  and  the  plume  becomes  the  dominant  source.  Also, 
if  the  aircraft  is  viewed  from  the  forward  hemisphere  or 
from  an  aspect  angle  at  which  the  tailpipes  are  not  visible, 
the  pl\ime  is  the  only  source  of  radiation  available.  The 
calculation  of  the  exact  radiant  intensity  of  the  plume  is 
extremely  difficult  since  both  temperature  and  emissivity 
vary  in  a  complex  manner  through  out  its  volume. 

The  exhaust  temperature  contours  for  a  typical  turbojet 
engine  with  and  without  afterburner  are  shown  in  Figure  4.1. 
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Figure  4.1  Exhaust  Temperature  Contours  for  a  JT4A  Turbojet 
Engine  with  and  without  Afterburner  (From  Ref  [1] ) 

We  can  see.  that  when  afterburner  is  turned  on,  the 
temperature  and  size  of  the  plume  increase  appreciably.  If 
we  integrate  over  the  entire  plume  in  order  to  estimate  its 
radiant  intensity,  the  radiant  intensity  will  be  several 
times  that  of  the  hot  tailpipe.  For  in-band  engineering 
calculations,  a  plume  can  be  approximated  as  a  graybody  with 
an  emissivity  of  0.9  [Ref  1]. 

The  steps  that  are  followed  to  calculate  the  radiant 
intensity  of  the  plume  model  is  shown  in  Figure  4.2. 
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Figure  4.2  Calculation  of  In-Band  Radiant  Intensity. 


B.  DEVELOPMENT  OF  THE  PLUME  MODEL 

The  assTomed  geometric  model  for  the  plume,  shown  in 
Figure  4.3,  is  an  ellipse.  This  is  used  to  model  a  three- 
dimensional  radiating  ellipsoidal  source.  The  required 
parameters  are  the  included  major  axis  the  minor  axis 

and  the  temperature  T„  that  is  assumed  constant 
throughout  the  plume.  The  area  of  the  ellipse  is 

Api=^-aprKi 


Figure  4.3  Plume's  Radiating  Source  Area. 

The  Stefan-Boltzmann  Law  is  again  used  to  predict  the 
radiant  emittance.  Here  it  is  applied  to  the  case  of  the 
plume: 

W{A,T)„=£oT„*  {IV -2) 

where  8  here  is  the  effective  emissivity  of  the  plume. 

Using  the  M-file  plankpl.m  for  the  particular 
temperature  with  the  specific  wavelengths  of  the  detector's 

band,  we  can  compute  the  total  emittance  that  lies 


40 


inside  the  detector's  band.  Next,  we  find  the  fraction  of 
the  emittance  that  lies  between  wavelengths  and  X^  as 


(lV-3) 


The  physical  model  for  the  pliame  source  is  assxamed  to  behave 
as  a  Lambertian  emitter  [Ref  2].  Following  the  standard  rule 
for  calculating  radiance  from  total  emittance  [Ref  2]  ,  we 
get 


AT 

ly  PI 


n 


(/V-4) 


The  total  radiant  intensity  is  then  given  by 


^  n  ~  ^pi '  ^Pi 


(IV -5) 


and  the  pliime  radiant  intensity  that  lies  inside  the 
detector's  band  is 


y  (a,  ,  —  J  PI  •  “Hpi 

(TV -6) 

For  later  reference,  the  composite  image  signal  associated 
with  the  plume  target  is  defined 

S^J(\,\)„  (tv-1) 


41 


42 
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INTEGRATION  OF  PYROPHORIC /PLUME  MODELS 


The  main  program,  recorded  in  Appendix  B,  consists  of 
two  parts.'  The  general  structure  of  the  main  program  is 
represented  in  Figure  5.1. 


Generate  S/N  Ratio  with  \ 
variation  in  parameters  I 
over  range  of  time  J 


j  Generate  images  at 
one  instant 


V 


Figure  5.1  General  structure  of  the  main  program. 


Once  the  inputs  have  been  specified,  a  selection  in  the 
main  program  permits  execution  of  either  the  first  part, 
which  is  the  generation  of  the  signal-to-noise  ratio  (S/N) , 
or  the  second  part,  which  refers  to  the  generation  of  the 
pyrophoric  flare  and  plume  images . 
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A.  CALCmATION  OF  S/N 

The  first  step  of  the  flowchart  in  Figure  5.2,  the 
conputation  of  the  time  radiant  intensity  distribution  of 
the  pyrophoric  burn,  was  covered  in  Chapter  III.  One  typical 
distribution  is  represented  in  Figure  3.4. 

In  the  second  step  of  the  flowchart,  the  in-band 
radiant  intensity  calculation  of  the  pyrophoric  material  for 
each  instant  of  time  t  is  performed.  This  is  found  by 
performing  a  multiplication  of  the  temporal  variation  in  the 
radiant  intensity  J(t),  given  by  the  equation  (III-9),  with 
the  fraction  of  the  emitted  power  ripj  that  lies  inside  the 
detector's  band  given  by  equation  (III-2)  . 

Next,  the  in-band  radiant  intensity  of  the  pyrophoric 
material  for  each  instant  of  time  t  is  obtained  from 

j{X^,X^Pf  =  j{t)  -TJpji  (^“l) 

and 


where  N  is  interpreted  as  a  composite  noise  interference  for 
the  missile  tracking  system. 
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The  in-band  radiant  intensity  calculation  of  the  plxome, 
as  given  by  equation  {IV-7)  in  Chapter  IV,  is  the  last 
component  needed  in  the  S/N  prediction. 

For  the  calculation  of  S/N,  a  situation  similar  to  that 
shown  in  Figure  5.3  is  assumed,  where  the  pliame  represents 
the  signal  and  the  pyrophoric  the  noise. 


{signal) 


(noise) 


Figure  5,3  S/N  Representation. 

Based  on  the  preceding  discussion,  by  taking  the  ratio  of 
the  result  of  equation  (IV-7),  which  represents  the 
composite  plxime  image  signal,  with  the  result  of  equation 
(V-2) ,  which  represents  the  composite  pyrophoric  image 
signal,  we  obtain  an  estimate  for  S/N: 


S  _  'l(^^^2)pi 


(I' -3) 
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After  the  estimation  of  S/N,  the  main  program  permits  a 
change  in  parameters  such  as  pyrophoric  temperature  T^j, 
plume  temperature  detector's  band  wavelengths  , 

pyrophoric  peak  time  tp,  or  radius  r^j^,  and  a  recalculation 
with  the  new  inputs  produces  a  new  value  for  S/N. 

The  second  part  in  the  main  program  refers  to  the 
generation  of  the  pyrophoric  and  plume  images.  The  main 
steps  of  the  procedure  which  are  required  to  create  the 
images  are  shown  in  Figure  5.4. 


Figure  5.4  Generate  Images  Structure. 
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B.  GENERATION  OF  IMAGES 

Pyrophoric  Image  Construction 

The  first  four  steps  of  the  flowchart  in  Figure  5.4 
have  been  already  created.  Specifically,  this  includes  the 
in-band  radiant  intensity  calculation  for  each  instant  of 
time  t.  This  calculation  is  done  in  the  first  part  of  the 
main  program  which  also  deals  with  the  S/N  generation.  Also 
covered  in  these  four  steps  is  the  conversion  from  the 
normalized  space  distribution  to  the  integrated  pyrophoric 
model  discussed  in  Chapter  III. 

To  construct  the  pyrophoric  image,  it  is  assumed  that  a 
situation  similar  to  Figure  5.5  is  applicable,  where  FOV 
represents  the  f ield-of-view  of  the  incoming  missile,  R  is 
the  distance  at  which  the  missile  will  detect  the  pyrophoric 
source,  and  the  nximber  of  pixels  that  have  to  be  used  in 

order  to  visualize  the  pyrophoric  image.  These  quantities 
are  some  of  the  inputs  which  are  specified  at  the  beginning 
of  the  main  program. 

From  this  assumption,  the  missile's  window  dimension 
is  calculated  as 

Rpidd  =^OV{rad)xR{m)  .  [m]  (’^^“4) 
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Figure  5.5  Pyrophoric  Image  Construction. 

To  specify  the  distance  of  the  pixels  from  the  center  of  the 
pyrophoric  image,  the  Pythagorean  relation  is  used 

Nr  =#17  (''-5: 

where  the  i,  j  indices  take  values  from  -(Np.^/2)  to  (N’pi^/2)  . 
Once  has  been  calculated  for  all  pixels,  the 

corresponding  radial  distances  are  predicted  from 


Next,  in  the  algorithm  an  index  is  generated  in  order  to 
assign  a  value  from  the  integrated-view  radiant  intensity, 
equation  (III-23)  ,  to  each  pixel.  This  index  is  deteonnined 
by 


1  = 


R 


pf 


[r^Ax  •max(r^o^)) 


(V-l) 


For  each  value  of  this  index  that  corresponds  to  the  inside 
of  the  integrated-view  radiant  intensity  array  index,  the 
corresponding  radiant  intensity  value  is  assigned 


(y_8) 

otherwise  the  zero  value  is  given  to  the  pixel.  Finally,  a 
normalization  is  performed.  The  value  of  each  pixel  is 
divided  by  the  summation  of  all  the  pixels  and  multiplied 
with  the  in-band  pyrophoric  radiant  intensity,  equation (V- 
1)  ,  The  normalization  rule  is  therefore  equivalent 
to 


C,, 


(y-9) 


where  the  time  dependence  is  implicit. 
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2.  Plume  image  construction 

From  the  flowchart  in  Figure  5.4,  the  only  necessary 
information  to  create  the  basic  plume  image  is  the 
calculation  of  in-band  radiant  intensity  of  plume  material, 
which  has  been  already  estimated.  Specifically,  the  in-band 
radiant  intensity  calculation  for  each  instant  of  time  t  is 
given  from  equation  (IV-7) . 

To  construct  the  plxome  image,  it  is  assumed  that  the 
geometric  conditions  represented  by  Figure  5 . 6  are 
applicable. 


Figure  5.6  Plume  Image  Construction. 
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Note  the  main  difference  between  Figure  5.5  and  Figure  5.6 
is  that  the  plume  is  represented  here  as  an  ellipse  instead 
of  as  a  circle.  The  FOV,  R,  and  are  the  same  quantities 

that  already  have  been  specified  for  the  pyrophoric  image 
construction. 

To  specify  the  location  of  the  pixels  in  the  plume 
image,  one  vector  r  is  created 

r  =  iA-e^  +  jA-  Cy  (V  -lO) 

where 

-  i,;j  denote  indices  that  take  values  from  -(N’p.^/2)  to 

(Npi*/2)  with  a  specific  step 

-  Cx,  Cy  denote  unit  vectors  along  the  major  axis  a  and 

minor  axis  b  ,  “ 

pi 


^  _  ^FieU 


(v-n) 


In  order  to  create  an  ellipsoidal  shape,  a  value  has  to  be 
assigned  to  the  pixels  that  belong  inside  the  ellipse's 
boundary 

(V-12) 


+^^<1 


2  '  L  2 


where  x,y  are  defined  as 


Taking  the  preceding  into  account  and  making  use  of 
equations  (V-11) ,  {V-12),  (V-13),  we  get 


(V-14) 


which  provides  a  convenient  rule  for  defining  the  interior 
part  of  the  plxime  in  the  computer  model . 

The  niomeric  values  which  are  assigned  to  the  pixels 

that  define  the  in-band  plume  radiant  intensity  are 

taken  from  pliime  model  and  equation  (IV-6)  .  After  the 
generation  of  the  ellipsoidal  shape,  a  normalization  takes 
place.  The  value  of  each  pixel  is  divided  with  the  summation 
of  all  the  pixels  and  nailtiplied  with  the  in-band  plume 
radiant  intensity  This  normalization  rule  is 

equivalent  to 


(V-15) 


The  similarity  of  equation  (V-15)  to  equation  (V-9)  for  the 
pyrophoric  image  is  apparent. 
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VI.  PRESENTATION  OF  S/N  CURVES 

As  mentioned  in  Chapter  V,  the  first  part  of  the  main 
program,  which  is  recorded  in  the  Appendix  B,  is  primarily 
concerned  with  the  generation  of  S/N.  Once  the  inputs  have 
been  specified,  the  value  of  S/N  can  be  extracted  for  each 
instant  of  time.  A  temporal  variation  in  S/N  during  the 
pyrophoric  burn  time  can  then  be  obtained.  Next,  by  changing 
parameters  such  as  pyrophoric  temperature  detector  band 

,  pyrophoric  peak  time  t^,  or  radius  the  program 

computes  S/N  curves  for  the  new  inputs. 

A  selection  of  these  curves  which  have  been  generated 
through  the  procedure  described  are  presented  in  this 
chapter.  A  graph  that  presents  an  overall  variation  of  S/N 
during  the  pyrophoric  burn  time  and  a  chart  that  records  the 
minimum  S/N  for  different  parameters  are  included  in  each 
figure . 

Specifically,  Figure  6.1  is  a  plot  of  S/N  as  a  function 
of  the  pyrophoric  function  time  for  different  values  of 
pyrophoric  temperature  T^j.  Results  were  generated  for 
Tp,=1000K,  tp=lsec,  and  r^=2m.  As  expected,  an  increase  in 
pyrophoric  temperature  results  in  a  decrease  in  S/N 

during  the  burn  period  of  the  pyrophoric  material.  This  is 
an  expected  result  because  it  is  known  from  the  Stefan- 
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Boltzmann  Law  that  the  emittance  increases  as  the  fourth 
power  of  the  temperature . 

Figure  6.2  is  a  plot  of  S/N  as  a  fxinction  of  the 
pyrophoric  function  time  for  different  values  of  radius  r 

HAX 

corresponding  to  the  instant  for  which  a  maximum  occurs 
in  the  radiant  intensity.  Results  were  generated  for 
'^pf“2000K,  Tpj=1000K,  and  tp=lsec.  We  observe  that  increasing 
the  radius  r^  significantly  decreases  S/N  between  the 
initial  and  final  time  of  the  pyrophoric  function.  This 
leads  to  the  conclusion  that  the  radiant  intensity  output  of 
the  pyrophoric  flare  increases  with  r  . 

HAX 

Figure  6.3  is  a  plot  of  S/N  as  a  function  of  the 
pyrophoric  burn  time  for  different  detector  bands.  Results 
were  generated  for  Tp,=2000K,  Tp,=1000K,  tp=0.8sec,  and  r^=2m. 
The  curve  with  the  lowest  values  corresponds  to  the  1-2 
micron  band.  This  band  exhibits  a  large  amount  of 
atmospheric  attenuation  and  is  not  a  practical  window  of 
detection  for  a  missile.  The  next  band  is  the  3-5  micron 
band,  which  has  less  atmospheric  attenuation.  This  band, 
S-Iong  with  the  8-12  micron  band,  exhibits  the  lowest  amount 
of  atmospheric  attenuation  and  are  practical  windows  of 
detection  for  a  missile.  Making  a  comparison  between  these 
two  bands,  we  see  from  Figure  6.3  that  the  curve  with  the 
lowest  values  corresponds  to  the  3-5  micron  band.  The  next 
lowest  band  is  the  7-9  micron  band.  Its  curve  has  lower 
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values  of  S/N  compared  with  the  8-12  micron  band,  but  again 
this  band  exhibits  a  large  amount  of  atmospheric  attenuation 
and  is  not  a  practical  window  of  detection  for  a  missile. 

Figure  6.4  is  a  plot  of  S/N  as  a  function  of  the 
pyrophoric  function  time  for  different  values  of  t^.  Results 
were  generated  for  Tpj=2000K,  Tp,=1000K,  and  r^=2m.  For  this 
case,  the  minimum  S/N  corresponding  to  each  are  the  same. 
The  focus  here  is  the  length  of  time  for  which  the 
pyrophoric  flare  will  achieve  intensity  that  can  capture  the 
missile's  seeker.  From  this  plot,  we  conclude  that  the  flare 
achieving  a  peak  radiant  intensity  at  a  later  point  in  time 
will  have  an  extended  windows  of  effectiveness.  This  point 
is  demonstrated  in  Figure  6.4  by  comparing  the  length  of 
time  S/N  is  below  0.5.  The  pyrophoric  flare  reaching  its 
peak  at  1.4  seconds  has  a  window  of  effectiveness  of  3,7 
seconds,  while  the  flare  with  tp=0.8  is  only  effective  for 
2.9  seconds.  The  former  is  "better"  by  approximately  25%. 
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S/N  Ratio 


S/N  Ratio 
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VXI.  PRESENTATION  OF  IMAGES 


The  second  part  of  the  main  program,  which  is  recorded 
in  Appendix  B,  consists  of  the  generation  of  the  pyrophoric 
and  plume  images.  The  procedure  which  is  followed  is 
described  in  Chapter  V. 

The  image  visualization  is  created  by  using  a  definable 
number  of  pixels.  The  command  "image"  in  MATLAB  creates  an 
image  graphics  object  by  interpreting  each  element  of  the 
pyrophoric  or  the  plume  matrix  which  has  been  created  as  an 
index  into  the  figure's  colormap.  Each  element  of  the  above 
matrix  specifies  the  color  of  a  rectilinear  patch  in  the 
image.  Along  with  the  command  "image",  the  command 
"CdataMapping/scaled"  is  used  to  scale  the  values. 

The  variation  in  radiant  intensity  for  the  pyrophoric 
image  can  be  described  better  with  two  colormaps.  The  "jet" 
colormap,  which  ranges  from  blue  to  red  and  passes  through 
the  colors  cyan,  yellow,  and  orange.  The  "gray"  colormap 
returns  a  linear  grayscale  colormap.  In  Figures  7.2  to  7.6, 
a  colorbar  is  used  to  show  the  current  color  scale. 

The  second  part  of  the  M-file  pyrof.m,  which  is 
recorded  in  Appendix  B,  is  used  to  provide  image  snapshots 
for  different  instances  of  time  for  the  pyrophoric  model  and 
for  the  plvome  model  which  is  constant  with  time.  For  example 
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see  Figure  7.1  which  was  generated  with  Tp3^=1000K.  A  sequence 
of  pyrophoric  images  on  a  "jet"  colormap  for  a  pyrophoric 
with  a  burn  function  of  seven  seconds,  burn  temperature  of 
Tpj=2000K,  detector  band  3-5  micron  band,  peak  radiant 
intensity  at  tp=1.5  seconds,  and  radius  at  tp  equal  to  two 
meters  are  shown  in  Figures  7.2  through  7.6  for  times  t=l, 
t=1.5,  t=3,  t=5,  and  t=7  seconds,  respectively.  Looking  at 

these  plots,  we  observe  that  the  ring  with  the  maximum 
^radiant  intensity  will  be  small  and  close  to  the  center  for 
small  times.  At  the  end  of  the  burn,  for  the  pyrophoric 
flare  the  circular  periphery  of  peak  radiant  intensity  has  a 
larger  size  but  a  lower  radiant  intensity  value. 

One  byproduct  of  the  time-space  modeling  of  the 
pyrophoric  burn  is  that  the  radiant  energy  profile  does  not 
remain  peaked  at  the  center.  This  point  was  previously  made 
in  Chapter  III  and  qualitatively  represented  with  the  time 
dependent  graph  of  Figure  3.8.  The  images  of  the  pyrophoric 
seen  in  Figures  7.2  to  7.6  demonstrate  this  property. 
This  is  a  byproduct  of  the  assumption  that  the  peak  in  the 
radiant  intensity  will  occur  for  tp>0  for  a  spherical  "shell" 
r=r^>0.  This  is  consistent  with  a  physical  model  of  an 
expanding  material  for  which  the  center  burns  first  and  then 
an  outward  radial  wave  of  combustion  follows. 
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Figure  7.1  The  Pliame  Image. 


■■■  ■ '■■C':-VjN.ufnbef  ofpixBis'.  ■'■••  \  ■■, '•■'V''' 

Figure  7.2  The  pyrophoric  image  at  t=l  sec. 
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Representation  of  the  images  may  give  the  impression  that 
setting  the  plume  in  the  center  of  the  pyrophoric  image  will 
allow  the  missile  to  detect  the  plume,  which  is  uncovered  as 
the  radial  wave  of  combustion  moves  outward.  This  is  not 
true  since  the  reticle  based  detection  system,  represented 
in  Figure  7.7,  converts  all  the  pixels  inside  the  FOV  to  one 
value . 


S{t)+N{t) 


Figure  7*7  Reticle  based  Detection  System. 


It  should  be  noted  that  the  S/N  predictions  in  Chapter 
VI  are  based  on  this  model.  In  particular,  using  a  reticle, 
the  missile  electro-optic  (EO)  process  converts  two- 
dimensional  (2D)  optical  image  to  a  one-dimensional  (ID) 
electronic  signal  [Ref  1] . 


68 


VIIZ. MULTIPLE  PYROPHORIC  FLARES 


The  main  computer  program  code  also  simulates  the 
dropping  of  multiple  pyrophoric  flares,  each  at  different 
instances  of  time  and  with  generally  different  parameters. 
These  parameters  include  pyrophoric  temperature,  T^j,  radius 
of  pyrophoric  at  the  peak  radiant  intensity,  r^y^^  and  the 
time  of  peak  radiant  intensity,  t^. 

Using  this  option,  a  scenario  can  be  generated  where  a 
pre-launched  unmanned  tactical  air  launched  decoy  (TALD) 
drops  a  sequence  of  flares  in  advance  of  a  manned  aircraft. 
The  missile  will  then  be  distracted  by  a  flare  generated  IR 
noise  curtain.  The  intent  is  to  provide  'cover'  to  the 
manned  aircraft  which  flies  into  a  safe-zone  corridor  a  few 
seconds  later.  This  is  represented  in  Figure  8.1. 

The  settings  on  the  different  parameters  of  the  ejected 
pyrophoric  can  be  chosen  according  to  the  conclusions 
reached  from  studying  Figures  6.1,  6.2,  and  6.4  in  Chapter 
VII.  In  particular,  an  increase  in  the  pyrophoric 
temperature  T^^  or  an  increase  in  the  radius  r^ 
significantly  decreases  S/N  between  the  initial  and  final 
time  of  the  pyrophoric  function.  Also,  a  flare  which 
achieves  peak  radiant  intensity  at  a  later  point  in  time 
will  have  an  extended  window  of  effectiveness. 
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The  radiant  intensity  for  four  pyrophoric  flares,  which  have 
a  function  burn  time  of  5s  and  are  launched  at  times  t=0s, 
t=0.5s,  t=ls,  and  t=1.5s  are  shown  in  Figure  8.2.  The  first 
flare  has  the  shortest  time  to  peak  radiant  intensity, 
tp=0.8s.  This  value  is  selected  in  order  to  capture  the 
missile  seeker  quickly.  Other  specifications  for  the  first 
flare  include  rjj^=2m  and  Tpj=2200K,  which  were  chosen  to  keep 
S/N  low.  The  second  and  third  flares,  which  are  launched  at 
t=0.5s  and  t=ls,  respectively,  have  a  longer  time  to  reach 
peak  radiant  intensity,  tp=1.2s,  in  order  to  create  an 
extended  window  of  effectiveness.  For  these  two  flares, 
Tpj=2000K  are  used.  These  have  been  selected  lower 
relative  to  r^  and  T^^  for  the  first  flare  since  there  is 
still  significant  radiant  energy  from  the  first  flare.  The 
last  flare  is  launched  at  t=1.5s  with  larger  values  of  r^ 
and  Tp£  relative  to  the  previous  two  flares,  r^=2m  and 
Tpj=2200K.  Higher  values  are  used  in  order  to  maintain  a  low 
value  of  S/N  by  compensating  for  the  reduction  in 
effectiveness  of  the  first  flare.  In  addition,  tp=1.2s  is 
used  in  order  to  have  an  extended  window  of  effectiveness. 

The  composite  radiant  intensity  of  the  four  flares 
versus  the  total  function  burn  time  is  shown  in  Figure  8.3. 
The  changes  in  slope  are  correlated  with  the  activation  of 
new  flares. 
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The  composite  S/N  verisus  the  total  function  burn  time 
is  represented  in  Ficfure  8.4.  Also,  in  this  plot  the  changes 
in  slope  are  correlated  with  the  activation  of  new  flares. 

The  S/N  curves  that  correspond  to  each  flare  are  shown 
in  Ficfure  8,5.  By  comparing  Figure  8.4  and  Figure  8.5,  it  is 
obvious  there  is  an  extension  of  the  length  of  time  S/N 
remains  below  some  threshold  (e.g.,  0.5)  with  the  use  of 
multiple  flares.  This  point  is  quantified  in  Table  8.1. 


Figure  8.1  Pre-la-unched  Unmanned  Tactical  Air  Launched  Decoy 
(TALD)  drops  a  Sequence  of  Flares  in  advance  of  a  Manned 
Aircraft . 
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Figure  8.2  The  Radiant  Intensity  Distribution  of  the  four 
Pyrophorics . 
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J(W/sr) 


time  (sec) 


Figure  8.3  The  Composite  Radiant  Intensity  Distribution  of 
the  four  Pyrophorics. 
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S/N  Ratio 


Fyophoric  Function  Time  (sec) 


Figure  8.4  The  Composite  S/N  versus  the  Total  Function  Burn 
Time. 
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time  (sec) 


Figure  8.5  S/N  ratio  Curves  for  Each  Flare. 


Table  8.1  Time  below  the  0.5  threshold. 


' ,v  i  f'.lare  number  . 

,  Time  below  th^  0.5  threshold  (sec) 

Flare  1 

~3 

Flare  2 

0 

Flare  3 

0 

Flare  4 

~3 

Composite  flare 

~5.5 
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In  order  to  demonstrate  how  these  results  apply  to 
tactical  decisions  that  need  to  be  made  in  the  field,  two 
examples  are  presented.  The  conditions  that  the  manned 
aircraft  be  given  safe  cover  while  being  viewed  by  a 
"stationary  missile"  can  be  stated  as 

/N  /  ^  ac 

where 

-At,s/„,  is  the  time  S/N  is  depressed  below  a  specified 
threshold 

"^dac  approximate  time  after  the  dropping  of  the  first 

flare  that  the  aircraft  enters  the  safe  corridor 

aircraft  velocity 


Equation  (VIII-1)  can  be  rearranged  as  follows; 


R, 


^dac 

/N 


iyiii-i) 


For  convenience  of  interpretation,  minimum  values  for 
(V,^/Rj)  are  plotted  in  Figure  8.6  versus  for  various 

values  in  the  time  delay,  t^^^. 
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Figure  8.6  Minimum  Velocity  Curves  for  various  Values  of 
Time  Delay  Parameter. 

Example  1 

Given  Vg^=186  m/sec  (=0.56  MACHl)  ,  'R=62  m,  t3^^=1.0s, 
what  is  the  minimum  allowed  time  to  keep  S/N  depressed  in 
order  to  provide  "safe"  cover? 

Step  1  Calculate  3s'^ 

Step  2  For  delay  t3^^=1.0s,  read  from  Figure  8.6,  1.5s. 

From  the  above  graph  and  reference  to  Table  8.1,  it  is 
obvious  that  any  one  of  the  flares  studied  would  be 
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sufficient.  However,  if  the  aircraft  delay  is  changed  to  3.0 
seconds  and  m/sec,  then  the  condition  on  minimum 

changes  to  about  4 . 1  seconds  and  more  than  one  flare  packet 
will  need  to  be  dropped. 

Example  2 

Given  At,s,„,=2s,  R,=62  m,  t^^=1.5s, 

what  is  the  minim-urn  velocity  of  the  aircraft  in  order  to  be 
provided  "safe"  cover? 

Step  1  For  At(g^^,=2s  and  t^^=1.5s,  read  from  Figure  8.6, 
2.0s-^ 

Step  2  Use  the  results  from  step  1  with  Rj=62  m  to  obtain 
{V„,„n,=124  m/sec  (=37%  MACHl) 

The  point  of  this  example  is  to  demonstrate  that  a  time 
limit  in  the  required  reduction  of  S/N  forces  a  predictable 
lower  limit  on  the  velocity  of  the  manned  aircraft. 

From  Figures  8.4  and  8.6  several  significant  and 
general  observations  can  be  deduced.  First,  because  the 
period  of  time  in  which  S/N  stays  below  a  given  threshold  is 
finite,  as  can  been  seen  from  Figure  8.4,  there  is  an  upper 
limit  on  the  delay  that  can  take  place  before  the  manned 
aircraft  enters  the  safe  corridor.  The  longest  this  delay 
can  be  is  the  total  time  interval  that  S/N  is  below  the 

threshold.  Note  from  Figure  8.6  that  each  curve  predicts  the 

minimum  required  velocity  of  aircraft  which  approaches 

infinity  as  the  elapsed  time  for  S/N  to  be  depressed 
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asymptotically  approaches  the  delay.  For  example,  the  curve 
with  a  delay  of  three  seconds  approaches  infinity  from  the 
right  as  the  elapsed  time  S/N  depressed  approaches  three 
seconds . 

These  curves  predict  in  general  that  a  greater  delay  in 
the  aircraft  entering  the  corridor  requires  a  greater 
velocity  of  the,  manned  aircraft  in  order  to  guarantee  safe 
passage.  Conversely,  the  curves  when  applied  in  reverse 
suggest  a  minimiim  required  time  for  the  depression  in  S/N. 
In  particular,  from  Figure  8.6  it  follows  that  the  greater 
the  delay  and/or  the  slower  the  velocity  of  the  manned 
aircraft,  the  longer  S/N  needs  to  stay  depressed  below  a 
specified  threshold. 
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IX.  CONCLUSIONS  AND  FUTURE  ENHANCEMENTS 


The  objective  of  this  research  was  to  model  the  impact 
of  dropping  pyrophoric  type  flares  which  are  in  the  field- 
of-view  of  a  prelaunch  missile.  It  is  assumed  that  these 
flares  are  launched  from  an  Tinmanned  tactical  air  launched 
decoy  (TALD) .  The  manned  aircraft,  which  needs  protective 
cover  and  will  pres'umably  benefit  from  the  flares,  will 
follow  into  the  corridor  a  few  seconds  after  the  start  of 
the  flare  sequence. 

In  order  to  reach  this  objective,  a  simulation  model  of 
the  time- space  radiant  intensity  distribution  for  a 
pyrophoric  flare  expendable  has  been  developed.  The 
performance  of  the  pyrophoric  flares  to  create  a  distraction 
for  the  missile  was  characterized  in  term  of  S/N  where  S  is 
the  radiant  power  of  the  plume  and  N  is  the  radiant  power  of 
one  or  more  flares.  Several  conclusions  from  this  study  are 
described  in  the  next  two  paragraphs . 

The  radiant  power  for  the  pyrophoric  flare  burn  in  the 
[3-5]  micron  band  was  foTind  to  be  significantly  higher  than 
what  occurs  in  the  [8-12]  micron  band.  It  was  also  found 
that  increasing  the  pyrophoric  tenperature  T^^  or  increasing 
the  radius  r^  significantly  decreases  S/N  between  the 
initial  and  final  time  of  the  pyrophoric  function.  Finally, 
a  flare  which  achieves  a  peak  radiant  intensity  later  in 
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time  has  an  extended  window  of  effectiveness;  i.e.,  a  longer 
window  in  time  for  which  S/N  remains  below  a  specified 
threshold. 

In  this  thesis  the  concept  of  dropping  of  multiple 
pyrophoric  flares  to  create  a  "safe"  corridor  for  manned 
aircraft  missions  was  evaluated.  In  particular,  multiple 
"packets",  lengthen  the  window  over  which  S/N  remains 
below  an  acceptable  threshold. 

Lastly,  it  was  noted  that  if  the  manned  aircraft  is  not 
moving  with  high  enough  velocity  there  will  be  a  problem  in 
generating  "safe"  cover.  From  this  concept,  a  rule  relating 
minimum  aircraft  velocity  to  both  the  delay  in  entering  the 
safe  corridor  and  the  time  S/N  is  depressed  was  presented 
and  illustrated  with  several  examples. 

A  nvimber  of  simplifying  assrimptions  were  applied  in  the 
study.  Future  work  on  the  simulation  model  could  focus  on 
one  or  more  of  these  assumptions.  For  example,  to  mention  a 
few  points,  atmospheric  absorption  was  neglected  in  this 
study.  Also,  the  pl\ime  was  ass\imed  to  be  a  constant 
temperature  graybody.  A  better  model  for  the  atmospheric 
absorption  would  be  to  consider  the  cold  CO^  absorption 
spectriim.  For  the  plxome,  consideration  of  the  hot  CO^ 
emission  spectrum  would  be  also  beneficial  (the  red  spike- 

®P^ke  effect)  .  Also,  an  experimental  data  based  model 
for  the  pyrophoric  burn  radiation  spectrum  could  be  used  to 
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refine  the  simulation  model.  Finally,  the  calculations  of 
S/N  could  be  based  on  more  complex  image  processing 
algorithms  that  are  indicative  of  modern  reticle  based 
missiles . 
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APPENDIX  A. 


COMPUTER  AND  ANALYSIS  VARIABLES 


.Anaiys-is' 

Brief  Description 

FOV 

•FOV 

Missile  f ield-of-view 

R 

R 

Range  between  missile  and  target 

Rfield 

^ield 

Range  inside  missile's  f ield-of- 
view  on  target's  plane 

Tpf 

T 

Pyrophoric  temperature 

Tpl 

Pliame  temperature 

lambda 1 

Initial  wavelength 

lambda 2 

Final  wavelength 

tl 

tl 

Initial  time  of  pyrophoric 
function 

t2 

tl 

Final  time  of  pyrophoric  fimction 

t_peak 

t. 

Pyrophoric  time  at  max  radiant 
intensity 

r_max_t 

r 

MAXt 

Radius  of  pyrophoric 
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ormax 

Radius  of  pyrophoric  at  max 
radiant  intensity 

rnorm 

•‘"NORM 

Normalization  radius  of 
pyrophoric 

Afl 

An 

Area  of  pyrophoric  at  max  radiant 
intensity 

Apl 

.  A,, 

Area  of  plume 

a 

Spi 

Major  axis  of  plume's  ellipse 
area 

b 

^Pi 

Minor  axis  of  pl\ame's  ellipse 
area 

vel 

Velocity  of  pyrophoric 

e 

e 

Average  emissivity 

s 

a 

Stefan-Boltzmann  constant 

Wf  lare 

Total  emittance  of  pyrophoric 

Wpliime 

W(X,T)^, 

Total  emittance  of  plume 

Wpfb 

^  '  ^2  )  pf 

Emittance  of  pyrophoric  in  the 
band  of  interest 

Wplb 

Emittance  of  pliime  in  the  band  of 
interest 

ratWfl 

^Pf 

Fraction  of  pyrophoric  emittance 
in  the  band  of  interest 
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ratWpl 

Fraction  of  pliame  emittance  in 
the  band  of  interest 

J_max 

Jkmc 

Peak  radiant  intensity  of 
pyrophoric 

J 

J(t) 

Radiant  intensity  of  pyrophoric 
function 

Jpfb 

Radiant  intensity  of  pyrophoric 
in  the  band  of  interest 

Jplb 

J  pi 

Radiant  intensity  of  plume  in  the 
band  of  interest 

fa 

f (t,a) 

Gamma  distribution 

G 

r(a) 

Gamma  function 

a 

a 

Parameter  (a)  of  Gamma 

Distribution 

fad 

f  (x;fl,a) 

Gaussian  space  distribution 

s 

a 

Standard  deviation  of  Gaussian 
distribution 

m 

JA 

Mean  of  Gaussian  distribution 

Q 

Q 

Amount  under  a  Gaussian 
distribution  curve 

scrJ 

3  (  ^NORM  '  ^  ^ 

Weighted  value  of  a  Gaussian 
distribution  in  normalized  space 
domain 

Intscr 

Jv(r„,t) 

Integrated  view  radiant  intensity 

Npix 

N'pix 

Number  of  pixels 

Cpf 

Pyrophoric  image  matrix 
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Cpl 

Ci, 

Plume  image  matrix 

Cpfn 

C ' 

Normalized  pyrophoric  image 
matrix 

Cpln 

Normalized  plume  image  matrix 

Ratio 

S/N 

S/N  ratio 
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APPENDIX  B. 


NATLAB  CODES 


This  appendix  contains  listing  of  all  MATLAB  input 
files  that  were  used  to  get  the  results  posted  in  Chapters 
VI, VII  and  VIII.  These  M-files  are: 

1.  Pyrof  .m, 

2.  Plankpf .m 

3.  Planlqpl.m 

4.  Dguadl.m 

5.  Dquadl.m,  and 

6.  Gamnaz.m 

Pyrof .m  is  the  main  computer  code  which  has  two  parts. 
The  first  part  generates  the  S/N  curves  for  one  or  multiple 
pyrophoric  with  different  parameters  and  the  second  part 
generates  the  pyrophoric  and  pl\ime  images. 


Planl^f.m  is  a  function  file  which  calculates  the 
Spectral  Radiant  Emittance  (SRE)  for  the  pyrophoric  at 
temperature  Tp^  (Kelvin)  . 


Plankpl.m  is  a  function 

file 

which 

calculates 

the 

Spectral  Radiant  Emittance 

(SRE) 

for 

the 

pliome 

at 

temperature  Tp^  (Kelvin)  . 

Dquadl.m  is  a  function 

file 

which 

calculates 

the 

Spectral  Radiant  Emittance 

(SRE) , 

for 

the 

plume 

at 

temperature  Tp^ (Kelvin),  inside  the  detectors  band. 
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Dquad2.in  is  a  function  file  which  calculates  the 
Spectral  Radiant  Emittance  (SRE) ,  for  the  pyrophoric  at 
temperature  (Kelvin),  inside  the  detectors  band. 

Gammaz.xn  is  a  function  file  which  calculates  the  GAMMA 
fiinction. 
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Program  1  Pyrof.m 


%%%%%%%%%%  MAIN  COMPUTER  CODE  %%%%%%%%%% 
%%%%Ref.  files: 

%  plankpf.m,  plankpl.m, 

%  dquadl . m ,  dquad2 ,  m ,  gamma  z . m 


format  long 

global  Tpl 
global  lambdal 
global  lambda2 
global  Tpf 
global  tl 
global  t2 

%%%%%  Inputs  %%%% 

program= input ( 'For  the  S/N  part  enter  1,  else  for  the  images  part  enter 
2 

%%%%%%%%%%%%%%%%%%%%%  First  Part 
if  program==l 

Tpl=input ( 'Set  the  Temperature  of  plume  (Kelvin):  '); 

Tp=input { 'Set  the  Temperature  of  Pyrophoric  (Kelvin):  '); 
lambdal=input  ( 'Set  the  value  of  the  lambdal  (\m)  :  '); 
lambda2 = input  { 'Set  the  value  of  the  lambda2  (lam)  :  '); 

tl=input ( 'Set  the  initial  value  of  the  time  for  the  Pyrophoric  (sec): 

' 

t2=input ( 'Set  the  final  value  of  the  time  for  the  Pyrophoric (sec) :  '); 

t_peakl = input  ( 'Set  the  value  of  the  t_peak  (sec):  '); 

rmaxp=  input  (' Set  the  value  of  the  rmax  (m)  :  '); 

t3=input ( 'Set  the  increment  step  between  tl  and  t2 :  '); 

tf l=input ( 'Set  the  dropped  times  of  pyroforic:  '); 

t=tl: t3 : t2; 
figure (1) 

for  z=l : length ( tf 1) ; 
tflare=tfl(z); 
if  length ( t_peakl ) ==1; 

t_peak- t_peakl ; 
else 

t_peak=t_peakl  ( z )  ; 

end 

if  length(Tp) ==1; 

Tpf=Tp; 

else 

Tpf=Tp  (z)  ; 

end 

if  length (rmaxp) ==1; 

rmax=rmaxp  ; 
else 

rmax=rmaxp  ( z )  ; 

end 

e=0.9;  %emissivity  of  Pyrophoric 

s=5 . 67*10^ (-12) ; %  Stef an-Boltzmann  constant  (cm) 
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Wpfb=dguad2  ( 'plankpf ' ,  lambdal,  lambda2)  ;%Einittance  between  the  two  lambda 
%(  W  cm^(-2)  ) 

Wf lare=s*Tpf^4;  %Stefan-Bolt2mann (total  emittance  over  all  wavelengths) 
Afl=4*pi* (rmax) ^2;  %  Area  Pyrophoric  at  t_peak 
sl=5 . 67*10^^  (“8 )  ;  %  Stefan-Boltzmann  constant  (m) 

J_max=  ( (e*sl*Tpf  ^^4)  /  (4*pi) )  *Af  1 ;  %  Total  Radiant  Intensity  (W/sr) 
ratWfl=  (Wpfb/Wflare)  ; 

m=  (tflare/t3) +1; 
a=t  _peak+l; 

G=gammaz(a);  %  Find  'G'  with  the  help  of  'gammaz'  function 
fa=  ( t . { t_peak)  .  *exp  (“t)  )  . /G;  %  gamma  equation 

J= (fa. *J_max) * (l/max(fa) ) ;  %  Normalized  to  the  total  Radiant  Intensity 

%J_max 

Js=J*ratWfl; 

[rl, r2] =size ( Js) ; 
cl=Js(l,2:r2) ; 
c= [zeros {l,m) ,cl] ; 
if  z==l; 

y=l: length (c) ; 
y=c; 

else 

y= [y, zeros (1, ( length (c) -length (y) ) ) ] +c; 

end 

tpl=0: t3 : ( ( length (c) -1) *t3) ; 

eval ( [ 'save  data' ,num2str (z) , ' .mat' ] ) ; 

eval ( [ 'load  data' ,num2str (z) , ' .mat' ] ) ; 

plot (tpl,c) 

grid  on 

xlabel ( ' time  ( sec ) ' ) 
ylabel ( ' J (W/sr) ' ) 
hold  on 
end 

figure (2) 

tcomp=0 : t3 : ( ( length (y) -1) *t3) ; 
plot ( tcomp,y) 
grid  on 

xlabel ( ' time  ( sec ) ' ) 
ylabel ( ' J (W/sr) ' ) 

%  Calculation  of  radiant  intensity  for  the  plume 

Wplb-dquadl  ( 'plankpl ' ,  lambdal ,  lambda2 )  ;  %Emittance  between  the  two  lambda 

%(  W  cm'"  (-2)  ) 

a=1.50;  %minor  axis  (m) 

b=3-50;  %major  axis  (m) 

Apl=pi*a*b;  % (Area) 
e=0.9;  %emissivity 

3=5.67*10^ (-12) ;  %Stefan-Boltzmann  constant 

Wplume=s*Tpl^4;  %Stefan-Boltzmann ( total  emittance  over  all  wavelengths) 

sl=5 .67*10'"  (-8)  ;%  Stefan-Boltzmann  constant  (m) 

ratWpl=Wplb/Wplume; 

Jpl= ( (e*sl*Tpl^4) /pi) *Apl*ratWpl;  %Radiant  Intensity  between  the  two 
%lainbda  (W/sr) 

Jplume=Jpl .  *ones  ( 1 ,  length  (y)  )  ; 

%  *********  S/N  Ratio  ******* 

Ratio= Jplume . /y ; 

figure (3) 

plot (tcomp, Ratio) 
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grid  on 

xlabel ( '  Pyrophoric  Function  Time  (sec)  ') 
ylabel('S/N  Ratio') 

figure (4) 

for  x=l : iength( tf 1) ; 
tf lare=tf 1 (x) ; 
if  length ( t_peakl ) ==1; 

t_peak==t_peakl ; 
else 

t  j)eak=t_peakl  (x)  ; 

end 

if  length (Tp) ==l; 

Tpf==Tp; 

else 

Tpf=Tp(x); 

end 

if  length(rmaxp) ==1; 

rmax==rmaxp  ; 
else 

rmax=rmaxp(x)  ; 

end 

%  Calculation  of  radiant  intensity  for  the  pyrophoric 
e=0.9;  %emissivity  of  Pyrophoric 
s=5.67*10'"(-12)  ;%  Stefan-Boltzmann  constant  (cm) 

Wpfb=dquad2  ( 'plankpf ' ,  lambdal,  lainbda2)  ;%Emittance  between  the  two  lambda 
%(  W  cm'^(-2)  ) 

Wflare=s*Tpf'^4;  %Stefan-Boltzmann (total  emittance  over  all  wavelengths) 
Afl=4*pi* (rmax) ^2;  %  Area  Pyrophoric  at  t_peak 
sl=5.67*10^(-8) ;%  Stefan-Boltzmann  constant  (m) 

J_max={  {e*sl*Tpf'"4)  /  (4*pi)  )  *Afl;  %  Total  Radiant  Intensity  (W/sr) 
ratWf  1=  (Wpfb/Wf  lare)  ; 

m=  ( tf lare/t3 ) +1 ; 
a=t_peak+l; 

G=gammaz(a);  %  Find  'G'  with  the  help  of  'gammaz'  function 
fa=(t.^(t_peak) .*exp(-t) ) ./G;%  gamma  equation 

J=(fa.*J_max)*(l/max(fa)) ;  %  Normalized  to  the  total  Radiant  Intensity 

%J_max 

Js=J*ratWfl; 

[rl,r2]=size(Js) ; 
cl=Js(l,2:r2) ; 
c2= [zeros (l,m) ,cl] ; 

Jplume=Jpl.*ones{l, length (c2) ) ; 
c=Jplume, /c2 ; 

tpl=0:t3: ( (length(c)-l)*t3)  ; 

eval ( [ 'save  datas' ,num2str (x) , ' .mat' ] )  ; 
eval ( [ ' load  datas ' , num2str (x) , ' .mat ' ] ) ; 
plot (tpl,c) 
grid  on 

xlabeK'time  (sec)  ') 
ylabeK'S/N  Ratio') 
hold  on 
end 

%%%%%%%%%%%%%%%%%%%%%  Second  Part 
else 

Npix=input ( 'Set  the  number  of  pixels  (it  must  be  an  integer  value  with 
half  an  odd  number) :'); 

tspf=input ( ' Set  the  time  (sec)  that  you  want  to  create  the  images:  '); 
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Tpl=input  ( 'Set  the  Temperature  of  plxome  (Kelvin):  '); 

Tpf=input { 'Set  the  Temperature  of  Pyrophoric  (Kelvin):  '); 

lambdal = input ( 'Set  the  value  of  the  lambdal  (urn):  '); 

lambda2=input  ( 'Set  the  value  of  the  lainbda2  (urn):  '); 

tl=input ( ' Set  the  initial  value  of  the  time  for  the  Pyrophoric  (sec): 

' )  ; 

t2=input ( 'Set  the  final  value  of  the  time  for  the  Pyrophoric (sec) :  '); 
t_^eak=input ( 'Set  the  value  of  the  t_peak  (sec):  '); 
rmax=input ( 'Set  the  value  of  the  rmax  (m) :  '); 
t3=input ( 'Set  the  increment  step  between  tl  and  t2 :  '); 

%  Calculation  of  radiant  intensity  for  the  plume 
Wplb=dquadl  ( 'plankpl ' ,  lambdal ,  lainbda2 )  ; 
a=1.5;  %minor  axis  (m) 
b=3.5;  %major  axis  (m) 

Apl=pi*a*b;  % (Area) 
e=0 . 9 ;  %emissivity 

s=5 . 67*10'^  (~12)  ;  %Stefan~Bolt2mann  constant 
Wplume=s*Tpl^4;  %total  emittance  over  all  wavelengths 
sl=5 , 67*10'"  (-8)  ;%  Stefan-Boltzmann  constant  (m) 

Jplb= ( (e*sl*Tpl^4) /pi) *Apl* (Wplb/Wplume) ;  %Radiant  Intensity  between  the 
%two  lambda (W/sr) 

FOV=( (pi*2)/180) ; 

Range=1.8*(10"3) ;  % (m) 

Rfield=FOV*Range;  % (m) 

%I  create  the  image  matrix  for  the  Plume 
Cpl=zeros (Npix+1 , Npix+1 ) ; 
for  il=-Npix/2 :l:Npix/2; 
for  jl=-Npix/2 :l:Npix/2; 
ipixl=il+( (Npix/2)+l) ; 
jpixl= j 1+ ( (Npix/2 ) +1 ) ; 

rl=((il)  .'"2./(a/Rfield)'^2)  +  ((jl)  .'"2  .  /  (b/Rfield)  ^2)  ; 
if  rl<=(Npix+l)'"2; 

CpKipixl,  jpixl)=Jplb;%  I  assign  the  value  Jplb  to  the  pixels 
%inside  the  matrix 
else 

Cpl (ipixl , jpixl) =0; %  I  assign  the  value  0  to  the  pixels  outside 
%the  matrix 

end 

end 

end 

suml=0; 

for  i=i: (Npix+1) ; 
for  j=l : (Npix+1) ; 

suml=  suml+  Cpl(i,j); 

end 

end 

newCpl=  Cpl/suml  *Jplb; 

Cpln=newCpl; 

%  Calculation  of  radiant  intensity  for  the  pyrophoric 

e=:0.9;  %emissivity  of  Pyrophoric 

s=5 . 67*10'"  (-12 );  %  Stefan-Boltzmann  constant  (cm) 

Wpfb=dquad2  ( 'plankpf ' ,  lambdal ,  lambda2 )  ;  %Emittance  between  the  two  lambda 
%(  W  cm'"  (-2)  ) 

Wf lare=s*Tpf ^4;  %total  emittance  over  all  wavelengths 
Afl=4*pi* (rmax) ^2;  %  Area  Pyrophoric  at  t_peak 
sl=5 . 67*10^ (-8) ; %  Stefan-Boltzmann  constant  (m) 
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J_inax=  {  (e*sl*Tpf  ^^4)  /  (4*pi)  )  *Af  1;  %  Total  Radiant  Intensity  (W/sr) 
ratWfl=  (Wpfb/Wf lare)  ; 

t=tl:0.1:t2; 
vel=rmax.  /t _peak; 
r_inax_  t = t .  * ve  1  ; 
rnorm=r_max_t  /rmax  ; 

%  I  create  the  time  distribution (gamma  distribution) for  a=t_peak+l 
a=t_peak+l; 

G=gammaz{a);  %  Find  'G'  with  the  help  of  'gammaz'  function 
fa=:{t.'^(t__peak)  .*exp(~t)  )  ./G;%  gamma  equation 

J= (fa. *J_max) * (l/max(fa) ) ;  %  Normalized  to  the  total  Radiant  Intensity 
% J_max 


%  I  create  the  normalized  space  distribution 
s=l;  %  standard  deviation 
for  tx=tl:t3:t2; 

fad=  (exp  (-'^(rnorm-  (tx/t_peak)  )  .  ^2  .  /  (2*s^2)  )  .  /  (sgrt  (2*pi)  *s)  )  ;  %Gaussian 
%distribution 
12=tx; 

fa2=(12'^  (t_peak)  *exp(-12)  )  /gammaz  ( t_peak+l)  ; 

J2=  (fa2*J__max)  *  (l/max(fa)  )  ;  %  Normalized  to  the  total  Radiant 
%Intensity  J_max 

Q2=(0.5) .*erfc( (-tx/ (t_peak*s) ) ./sqrt (2) ) ;  %  Integral  from  zero  to 
%infinity  for  a  Gaussian  distribution  with  mean  at  tx/t_peak 
scrJl=fad*J2/Q2;  %The  weighted  values  for  every  Gaussian  distribution 
%with  mean  at  tx/t_peak 
end 

for  z=l : length ( tspf ) ; 
ts=tspf (z) ; 

fay=  ( exp  ( -  ( rnorm-  (  ( ts/t__peak)  )  )  .  ^2  .  /  (2*s^2)  )  .  /  (sqrt  (2*pi)  *s)  )  ; 
ls=ts; 

fas= (Is^ (t_peak) *exp(-ls) ) /gammaz (t_peak+l) ; 

Js=  (fas*J_max)  *  (1/max (fa) )  ; 

Qs=(0.5) .*erfc( (-ts/ (t_peak*s) ) ./sqrt (2) ) ; 
scrJ3=fay*Js/Qs; 

%I  create  the  integrated  Pyrophoric  model 
Np=30; 

[rl,ci]=size(rnoim)  ; 
romax=max  ( rnorm)  ; 

Nro=cl-l; 

Intscr2=l :Nro+l; 

ro=l:Nro+l;  %  ro  is  the  distance  from  the  center  of  projection 
Dro=romax/Nro  ; 
for  J=l:Nro+l; 
sum=0 ; 

ro(J)=(J-l) ,*Dro;  %  for  every  J  pick  up  one  ro 
xmax=sqrt( (romax) .^2-(ro(J) ) .^2) ; 

Dx=xmax . /Np ; 
for  n=l:Np+l; 
x= (n-1) . *Dx; 

r=sqrt  (x.'"2+(ro(J)  )  .^2)  ; 
kl=  (r.  /max  (rnorm)  )  *cl; 

k=round(kl);  %  index  k  goes  to  closest  integer 
if  k>cl 
a=0; 

sum=sum+a; 
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elseif  k==0; 
a=0; 

s\ani=sum+a  ; 
else 

siim=sum+scrJ3  (k)  ; 
end 

end 

Intscr2  (J)  =siain; 

end 

%I  create  the  image  matrix  for  the  Pyrophoric  model 
FOV=( (pi*2)/180)  ; 

Range=1.8*(10''3)  ;  %  (m) 

Rfield=FOV*Range;  % (m) 

[rl,cl] =size<rnorm) ;  %  cl  takes  the  value  of  size (r_max_t) 
Cpf=zeros(Npix+l,Npix+l) ; 
for  i=-Npix/2 : 1 :Npix/2 ; 
for  j=-Npix/2 :l:Npix/2; 
ipix=i+( (Npix/2)+l) ; 
jpix=j+{ (Npix/2)+l) ; 

Nr=sqrt(i.''2+j  .■'2)  ; 

R= (Rf ield/ (Npix+1) ) *Nr; 

Ll=  (R/  (max(morm)  *rmax) )  *cl; 

L=round(Ll) ; 
if  L>cl 

Cpf ( ipix, jpix) =0 ; 
elseif  L==0; 

Cpf ( ipix, jpix) =0 ; 
else 

Cpf { ipix, jpix) =Intscr2 (L) ; 

end 

end 

end 

Jts=Js*ratWfl; 
siim=0 ; 

for  i=l: (Npix+1) ; 
for  j=l : (Npix+l) ; 

sim=  sum+  Cpf  ( i ,  j )  ; 

end 

end 

newCpf=  Cpf /sum  *Jts; 

Cpf n=newCpf ; 

%  *********  s/N  Ratio  ******* 

Ratio=Jplb/Jts 

eval ( [ 'save  data' ,num2str (z) , ' .mat' ] ) ; 
end 

for  z=l : length ( tspf ) ; 

eval ( [ 'load  data' ,num2str (z) , ' .mat' ] ) ; 
figure (z) 

image (Cpfn, 'CDataMapping' , 'scaled') 
end 

figure (z+1) 

image(Cpln, 'CDataMapping' , 'scaled' ) 
end 
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Program  2 


Plankpf.m 


%  This  function  file  calculates  the  Spectral  Radiant  Emittance  (SRE) 
%  for  the  pyrophoric  at  temperature  Tpf (Kelvin). 

function  Y  =  plankpf (Lambda, T) 

global  Tpf 

Cl=3.7415e4; 

c2=1.4388e4; 

y  =  cl  ./(Lambda  .^5  .*  (exp(c2  ./(Lambda  .*  Tpf) )-!)); 

%  Y  will  have  units  of  watts/ (cm'^ 2 -micron) 


Program  3  Plankpl.m 

%  This  function  file  calculates  the  Spectral  Radiant  Emittance  (SRE) 
%  for  the  plume  at  temperature  Tpl (Kelvin) . 

function  Y  =  plankpl (Lambda, T) 

global  Tpl 

cl=3.7415e4; 

C2=1.4388e4; 

Y  =  cl  ./(Lambda  .^"5  .*  (exp(c2  ./(Lambda  .*  Tpl)  )-!)); 

%  y  will  have  units  of  watts/ (cm'^2 -micron) 


Program  4  Dquadl.m 

%  DQUADl  Numerically  integrates  an  expression  using  QUADS. 

function  [ I ] =dquadl ( Fun , al , a2 , a3 , a4 ) 

var  =  ' X ' ; 
tol=le-3; 
if  nargin<3, 
help  dquadl 
return 
end 

if  nargin  ==  3, 
xmin=al ; 
xmax=a2 ; 

elseif  nargin  ==  4, 
if  (isstr (al) ) , 
var=al; 
xmin=a2 ; 
xmax=a3 ; 
else 
xmin=al ; 
xmax=a2 ; 
tol=a3; 
end 
else 

var=al ; 
xmin=a2 ; 
xmax=a3 ; 
tol=a4; 

end 

I=quad8  ( 'plankpl'  ,xmin,xmax,  tol)  ; 
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Program  5  Dquad2.m 

%  DQUAD2  Niimerically  integrates  an  expression  using  QUADS . 

f unc t ion  [ I ] =dquad2 ( Fun , al , a2 , a3 , a4 ) 

var  =  'X'; 

tol=le-3 ; 

if  nargin<3, 
help  dquad2 
return 

end 

if  nargin  ==  3, 
xmin=al ; 
xmax=a2 ; 

elseif  nargin  ==  4, 
if  (isstr (al) ) , 
var=al ; 
xmin=a2 ; 
xmax=a3 ; 
else 
xmin=al ; 
xmax=a2 ; 
tol=a3; 
end 

else 

var=al ; 
xmin=a2 ; 
xmax=a3 ; 
tol=a4; 

end 

I=guad8  ( 'plankpf '  ,xinin,xmax,  tol)  ; 


Program  6  Gammaz.m 

%  GAMMAZ  gamma  function  that  also  allows  for  complex  arguments,  unlike 
%  MATLAB's  GAMMA  function.  For  real  arguments,  GAMMAZ  gives  the 

%  same  results  as  GAMMA  (within  numerical  error) . 

function  [f]  =  gammaz(z); 

[m  n]  =  size (z)  ; 
f  =  zeros  (m,n); 
cof  =  [76.18009172947146 
-86.50532032941677 
24.01409824083091 
-1.231739572450155 
0.1208650973866179e-2 

-0 . 5395239384953e-5] ;  %  6  coefficients  in  series  expansion 

for  icol  =  l;n,  %  do  one  column  at  a  time 

zz  =  z(: ,icol) ; 

zp  =  ones { 6 , 1) *zz '  +  [1 : length{ cof )]' *ones (1, length (zz) ) ;  %  vectorize 
ser  =  (cof '*(l./zp)+l. 000000000190015)  '  ; 
tmp  =  2Z+5.5  -  (2Z+.5) .*log(z2+5.5) ; 

Ingamma  =  -tmp  +  log(2.5066282746310005*ser./2z); 
f(:,icol)  =  exp ( Ingamma ) ; 
end 
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APPENDIX  C. 


MATHEMATICAL  DERIVATIONS 


1.  Solution  for  t^sa-l 

The  pdf  of  a  standard  Gamma  RV  is  given  by  equation  (III- 
7)  in  Chapter  III.  For  convenience  it  is  reproduced  below 

/(';«)  =  (C-1) 

f>0 
a  >0 

Now  by  setting 

with  r(a)  =  (a-l)!9i=0  for  a>0 
we  get 

r(fl) 

^[(a-i)f“"V'+r'^'*{-l)e"'J=0=>a  =  (+l 


(C-2) 


Now,  for  t=tp 

tp=a-l 


(C-4) 
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2.  Positive  side  integration  of  Gaussian  distribution 

The  pdf  of  a  Gaussian  RV  is  given  by  equation  (III-IO)  in 
Chapter  III,  For  convenience  it  is  reproduced  below 

1  Jx-M)/  ^  -oo<X<oo 

f(x,^,(T)  =  -j=—-e  2.^  -oo<^<oo  mean  (CS) 

0<(T  ^  tan  dard  deviation 


The  above  equation  for  ^.=0  and  a=l  becomes 


/(r;0,l)=  .  ^  -e  ^  -o«><jc<oo 
\2-7r 


(C-6) 


To  find  the  area  under  the  curve  from  x  to  infinity  the 
following  integral  must  be  evaluated: 

In  terms  of  the  complementary  error  function 

(C-8) 

Now,  if  the  continuous  random  variable  X  has  mean  and 


(C-7) 


Standard  deviation  a  then 


(C-9) 

(7 

is  a  Gaussian  random  variable  with  ji=0  and  0=1,  and  the  area 
under  the  curve  from  z  to  infinity  is 

e(z;0,l)  =  -^r 

From  equation  {C-8),  we  have 

Qiz)  =  )^'erfc  (C-ll) 

The  area  under  the  curve  from  zero  to  infinity  is  obtained 
by  setting  x=0  in  the  previous  equation  to  obtain  [5] 


(C-12) 
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